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Abstract. We present numerical results for pure SU(2) Yang-Mills theory in four space-time dimensions 
using a novel algorithm based on dually transformed variables. The simulation makes use of a recently 
derived 0(j' 4 ) algorithm for the dual vertex amplitude and a dual Metropolis algorithm that generalizes the 
one recently developed for three dimensions. The dual algorithm is validated against the equivalent model 
using conventional variables over a range of couplings, spin cut-offs, and lattice sizes. We consider a lattice 
size up to 8 4 , where the problem of negative amplitudes renders the simulation results excessively noisy even 
at a relatively low /3 (starting at about j3 = 1.8). In conclusion, we survey some approaches to addressing 
the sign problem and increasing the efficiency of dual computations. 

1. Introduction 

Exact duality transformations for statistical mechanical systems with many degrees of freedom have 
provided a rich source of insights and alternative computational methods over the years [31[T7] . As lattice 
regularized quantum theories share much of the essential structure of lattice statistical mechanical models, 
it is natural to investigate analogous duality transformations. Indeed, applications of duality to quantum 
fields on the lattice continues to grow; recent examples include [51 I29H31] . In the present work, we shall focus 
on pure Yang-Mills theory as we expect it will be instructive to construct dual algorithms in this case before 
proceeding to dynamical fermions in the dual. 

For lattice gauge theory (LGT) with Abelian groups such as U(l) and Z2, it is well understood how to 
construct dual models useful for practical computations 22-24,26 . However, dual simulations with non- 
abelian gauge groups have presented a greater challenge. The main difficulties have been the definition of the 
dual model in terms of practically computable amplitudes and the construction of an ergodic set of moves 
that can be used with a Markov chain Monte Carlo method (e.g. the Metropolis algorithm). 

Following some initial results on amplitudes for non-abelian lattice models in the specific case of D = 3 
and G = SU(2) in [flU], a procedure for obtaining non-abelian dual models in four and higher dimensions 
on the lattice was given by Halliday and Suranyi [15) . More recently, a new mathematical framework for 
understanding the non-abelian dual transformation was constructed, the lattice spin foam formulation |20[ 
I21j . With regard to numerical simulations, the spin foam approach is helpful for computing amplitudes 
(particularly in the four-dimensional case [8]) as it gives an explicit prescription for the dual amplitudes in 
terms of spin networks. These spin networks can in turn be reduced in complexity from their initial definition 
by the use of diagrammatic recoupling methods [4UTTJ], thus avoiding the typically explosive proliferation of 
tensorial algebra that would otherwise ensue. The spin foam construction also makes more transparent the 
distinction between irrep labels associated to plaquettes and intertwiner labels that are naturally associated 
to edges, as well as the admissibility constraints between these quantities necessary for a non-zero amplitude. 
In [7], these features were applied to the three dimensional case with G = SU(2), resulting in an explicit 
algorithm that was validated against conventional results. Further details on the development and current 
status of this area can be found in [9] . 

The present work may be viewed as a follow-up in four dimensions to the previous results in three 
dimensions that includes [7] as well as earlier work by Hari Dass et al. [T2lI14| . In the three-dimensional 
case, an explicit amplitude was known from [HE] and the primary challenge was finding an ergodic set of 
moves. In the present case, the four-dimensional (hyper-cubic) vertex amplitude is a non-trivial function of 
48 spin variables, and it wasn't until recently that a practical algorithm for this vertex amplitude [5] was 
available. 

In this paper, we construct a set of ergodic moves (detailed in the Appendix) which when combined with 
the amplitude of [8] allow us to perform the first Metropolis algorithm simulations of dual SU{2) in four 



dimensions. We emphasize that the primary purpose of the present work was to validate the dual algorithm 
and its implementation against a conventional algorithm, in order to understand its performance over a 
range of couplings and to evaluate the severity of the sign problem in a simple setting. Efforts to extend the 
methods to more observables such as Wilson or Polyakov loops are presently underway by the author. 

This paper is organized as follows. In Section 2, we briefly summarize the form of the dual model for 
an SU(2) lattice gauge theory and define the observable to be computed. In Section 3, we describe how 
a conventional lattice code was constructed and used to compute expectation values of an appropriately 
chosen effective observable. Section 4 summarizes the results obtained with the dual algorithm. The paper 
is concluded with some discussion of the sign problem and other features of the dual algorithm in Section 5. 
Details of the moves that define the dual algorithm are provided in the Appendix. 

2. The Spin Foam Simulations 

In the present section we shall briefly review some essential features of the dual spin foam model in order 
to describe the simulations performed and results obtained. The reader is referred to [7] for a derivation of 
spin foam duality motivated by simulations with G = SU(2) in three dimensions; its adaption to the present 
case of four dimensions can be found in [S]. A general derivation of spin foam models can be found in [21J. 

In a dual spin foam model, the continuous group- valued variables of the conventional theory are replaced 
by discrete labels. The plaquettes of the lattice are labelled by irreducible unitary representations of SU(2), 
typically taken as the half-integers (although whole integers are sometimes used in diagrammatic methods). 
For a given plaquette labelling, each edge carries an intertwiner label, corresponding to one of a complete 
basis of maps that intertwine the irreps on plaquettes that intersect the edge. If an appropriate set of 
conditions is satisfied, than a particular labelling of plaquettes by irreps and edges by intertwiner labels is 
said to form an admissible spin foam. Further details on these variables and the admissibility constraints 
among them are given in the Appendix. 

As our present motivation is to evaluate the correctness and efficiency of the dual algorithm relative to the 
conventional formulation, we chose an observable that is particularly simple to compute in the dual model. 
Such an observable is provided by the average value of the spin at a plaquette p* : 

[) {Jr)D ~ E /e ^(/) ' 

where JeF are the admissible spin foams, and A(f) is the dual amplitude, as reviewed in Appendix |A"1 
The observable j p is the spin label assigned by spin foam / to plaquette p. The subscript D on the LHS 
indicates evaluation of the observable against the dual (as opposed to conventional) amplitude A(f). In 
Section [3l we will describe how to define and compute the same observable within the conventional LGT 
framework. 

For concreteness, we define here precisely which amplitudes were used to generate the results given in 
Section [U In terms of the spin foam amplitude A(f) defined in equation © of the Appendix, the vertex 
and edge amplitudes A v and A e are those defined in [S]. The plaquette amplitude A p used comes from the 
heat kernel action, defined as follows: 

Definition 2.1 (Heat Kernel Action). Let G = SU(N) and g 6 G. Let f(g,t) be the solution to the 
differential equation 

(2) A/ M = ^M 

at t = with initial condition f(g, 0) = S(g). Then the action S(g) defined implicitly by 
(3) 

is known as the heat kernel action [19l. 



S(g) 



f(g,t) 

f(I,t) 



It can be shown that the character expansion of the heat kernel action for SU(2) is of the form 
(4) e- p8 to = YX2j + l)e-%xj(9) 
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where Cj is the quadratic Casimir eigenvalue associated to the jth representation; for SU(2), Cj — j(j + 1). 



As the xid) are absorbed into the vertex and edge amplitudes of the dual model, by inspection of (U|) the 

plaquette amplitude A p in equation (J8J) associated with the heat kernel action is thus A p (j p ) = (2j p + l)e~^ s &~ ; 
the reader is referred to [7] for further details. 

3. The Conventional Simulation 

As we are working with a mathematically exact transformation between the conventional and dual forms 
of the lattice gauge theory, it is possible to validate the correctness of the dual algorithm by computing the 
same expectation value on both sides of the duality. We believe it is important to do so for several reasons. 
The dual vertex amplitude we make use of is relatively complex in both its derivation and hnal form [5]. 
Assuming that a dual amplitude has been correctly derived and implemented, one would also like to verify 
that the algorithm constructed (see the Appendix) is in fact ergodic and capable of realizing convergence of 
expectation values with comparable time and resources to that of a conventional code. Agreement within 
statistical error of conventional and dual results, while not constituting proof of a correct dual algorithm, 
strongly indicates its correctness. Such testing is particularly valuable as the dual code is further optimized. 

Given an observable and its expectation value defined in terms of the dual model, one can derive the 
effective observable whose expectation value is equal when evaluated with respect to the original conventional 
model. Explicitly, we seek a function Oj p (g e ) such that 

(5) (J P ') D = {O ip ,) c =\ []Jdg e O jp ,( ge )e^ s ^ 

J eeE 
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where we have introduced a character expansion in the second line to better show the relation with the dual 
observable. We next recall from Section [2] the dual observable of interest 

(6) (jp*) D = % Yl V II M{J}v,{i}v) II M{j}e,ie) II MJp) 
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Based on the above we introduce the following proposal for Oj p (g e ): 

( £ 3> eZ+ Jp' dim(j p . ) CjXjp . (ffp» ) \ 

The above can be checked by substituting into the second line of equation ([5]) and comparing to the second 
line of equation ([6]). 

Given the above, it is useful to be able to readily compute the character expansion components Cj of the 
action in order to compute the "re- weighted" numerator of ([7]). The heat kernel action used in the present 
simulations has a simple form (see Definition 12. ip in which character coefficients can be rapidly evaluated. 
We note in passing that the division by part of the original amplitude as is done in ([JJ in general may produce 
numerical problems, as convergence is better when observables are well correlated with the amplitude. In 
practice, for this particular case we have found good convergence (low variation of independent runs) of 
simulations for this effective observable. 



4. Description of Simulations and Results 

Several sets of computations were performed in order to confirm the correctness of the algorithm and eval- 
uate the computational performance of the dual simulations as currently implemented. These are described 
in the subsections below; we first make some general comments. 

Prior to computing the observable (j), the conventional code was validated against published curves 
appearing in [18] , which we selected because those simulations used the same heat kernel action employed 
in the present work. 
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Figure 1 . Comparison of conventional with spin foam code for 2 4 lattice at j c = 1 



Due to the computational cost of the vertex amplitude at higher spins, for the simulations described below 
the series appearing in (|5|) are cut off at a relatively low spin j c ; however, because the cut-off is applied to 
both conventional (as a truncation of the character expansion) and dual forms, a direct comparison can still 
be made, and posible errors in the (considerably more complex) dual code can be detected. While low cut-off 
simulations can be quite accurate in the strong coupling regime, they begin to break down at the weaker 
coupling scales of interest. Further optimizations of the vertex amplitude underway by the author should 
allow for higher spin cut-offs. 

We now turn to a description of the numerical simulations performed on both sides of the duality that 
establish the correctness (within statistical error) of the dual algorithm. Error bars of 3cr are shown, where 
a is the standard deviation based on the average and variance of over the different simulation runs. 

4.1. Small lattice, small cut-off case. Our first simulations were done on four-dimensional hypercubic 
lattice of 2 4 vertices and a cutoff of j c — 1 (in units of half-integer spin). The very small lattice size and low 
cut-off allowed relatively rapid simulations on both sides of the duality in order to establish agreement. 

For the conventional simulation, updating was done in "sweeps" in which each edge of the lattice was 
updated in sequence. A total of 10 7 sweeps were done at each value of (3. For the dual simulations, a total 
of 2 x 10 9 moves were applied; these moves are described in the Appendix. 

4.2. Small lattice, cut-off sensitivity analysis. Further comparisons of dual and conventional simula- 
tions were made at increased spin cut-offs of j c = | and j c = 2. At these cut-offs, 1 x 10 9 samples were used 
in the dual simulations on the minimal (2 4 ) lattice. The results are shown in Figure [5] and again confirm the 
correctness of the dual algorithm to within statistical error of the equivalent conventional algorithm. Error 
bars at weak coupling are larger at weak coupling, likely due to the smaller number of samples relative to 
the lower variance results obtained for j c = 1. 

4.3. Large lattice, small cut-off. The dual algorithm was also investigated on larger lattice having di- 
mensions 8 4 at a spin cut-off of j c = 1. At stronger coupling ((3 — 1.8 and lower) the dual simulations have 
high quality estimates with little relative error. However, for (3 and greater than 1.8 considerably larger 
errors begin to appear, rendering the results essentially unusable. As expected, one clear source of error is 
in the expectation value of the sign, which at higher beta tends to smaller values. The error at (3 = 1.8 on 
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Figure 2. Comparison of conventional with spin foam code for 2 4 lattice at j c = | 

the 8 4 lattice was 22% and the sign expectation was of order 0.01. For comparison, on the 2 4 lattice the 
error at a comparable (3 was 3% and the sign expectation value was of the order 0.1. 

5. Conclusions 

We have described the construction of a dual algorithm for SU{2) Yang-Mills in D = 4, based on a set 
of ergodic moves (suitable for use with a Metropolis or other algorithm) defined in the Appendix, and an 
efficient form for the vertex amplitude, derived in [8]. 

Using a 2 4 lattice as a simple test case, this algorithm was validated for a range of (3 at spin cut-offs 
of j c = 1 and j c — |. In this low-spin, small-lattice, intermediate coupling regime our algorithm exhibits 
statistical agreement with conventionally obtained results. While this is clearly an impractically small size, it 
represents a milestone as the first dual simulation of its kind in four space-time dimensions for a non-abelian 
gauge theory. 

Towards strong coupling (/? of order unity), accurate results are possible even with the low cut-offs used 
(as large spins are penalized heavily by the amplitude in this regime). However towards weak coupling 
an increased spin cut-off is required as configurations with larger spin begin to dominate. Because the 
computational cost of the dual vertex amplitude used increases with higher spins, raising the cut-off has the 
effect of slowing the rate at which samples are accumulated. Thus, a way of more efficiently obtaining high 
accuracy results at weaker coupling would be to speed up the vertex amplitude at high j, either through an 
improved algorithm, improved implementation, or some combination of both. Note that, as demonstrated 
in the D = 3 case [7] , different choices of intertwiner bases and reductions of the spin network via recoupling 
moves lead to algorithms of different performance. Therefore exploring these different choices is a natural 
way to seek out more efficient vertex amplitude algorithms. Intelligent caching of often used intermediate 
calculations and the use of recursion relations that exist for 6j symbols are also being investigated by the 
author. 

A serious challenge for the new algorithm was encountered in increasing the lattice to 8 4 in size. For this 
lattice size the sign expectation value had a considerably more rapid fall-off with increased 0. As the overall 
sign of amplitude is the product of sign factors local vertex amplitudes, a decrease in the sign expectation 
value for any given [3 can be generically expected as the lattice size becomes large, as was discussed in by 
Dass et al. in [14] . The rapid increased in error on data points for (3 = 1.8 and greater exhibits has allowed 
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us to measure for the first time the extent of the sign problem for the four-dimensional case (within our 
current approach). 

A number of approaches are being pursued by the author to mitigate the problem of fluctuating signs. One 
general approach would be to replace the 6j symbols that appear in the amplitude by their approximation 
(valid at large spins) in the form of well-known the Ponzano-Regge asymptotic formula [35] • While the 
Ponzano-Regge approximation itself contains fluctuating signs, because the expressions are in closed form 
they may useful as predictors of sign for the original 6j symbols and specifically to form a stationary phase 
approximations, in which configurations of stationary phase are used rather than the full ensemble [llj . A 
possible problem with this approach is that the regime where asymptotic formulae for spins may correspond 
to unrealistically large j3. 

Another approach to sign problems (that could work directly with the original vertex amplitude) is to 
identify a sector of the original theory in which some form of cancellation occurs in a controlled fashion 
(leading to a positive or at least better behaved sector for simulation). This has been achieved for certain 
types of fermion models, under the rubric of meron cluster and related methods; in some cases, configurations 
with exactly equal weight and opposite sign can be identified and removed a priori from the simulation and 
the remaining ensemble sampled in an efficient manner [BJ. 

A general lesson for the sign problem is that while general methods that work in all instances don't exist 
in tractable form [28], by using details of how differently signed configuration occur, it is possible to alleviate 
and in some casees completely remove the problem in a particular theory. A final point to make is that a 
useful solution does not necessarily have to work for arbitrarily high /?, but rather mitigate the sign problem 
to the point of being practical for j3 used to study length scales of particular interest. 

In summary, the present work has demonstrated for the first time that dual simulations of Yang-Mills 
in four space-time dimensions are feasible and provides a concrete, tested algorithm for carrying them out. 
The main barriers to extending this approach to large lattice sizes and weaker coupling are the emergence of 
a sign problem (the severity of which depends on lattice size) and the computational expense of evaluating 
the vertex amplitude; some possible approaches to addressing these were identified. 

Acknowledgement. The author would like to thank Dan Christensen for valuable discussion influencing this 
work. This work was made possible by the facilities of the Shared Hierarchical Academic Research Computing 
Network (SHARCNET) 

A. Dual algorithm for G — SU(2) in four dimensions 

As in the conventional framework of LGT, the setting for the present algorithm is a hypercubic lattice 
with periodic boundary conditions such that the lattice has the topology of a 4-torus (T 4 ). The partition 
function for the theory can be expressed as a sum over dual amplitudes as follows: 

(8) z YM = e A (f) = e n MUh, {i}v) n May* q n mj p i 

feF„ feF v£V e£E pGP 

where Fq is the set of lattice spin foam configurations (discussed further below). Each configuration labels 
a plaquette or edge by a quantum number. As can be seen from equation (|8|). the amplitude assigned to 
each spin foam factors into a product of contributions from local amplitudes associated to vertices, edges, 
and plaquettes of the lattice. Recently, an explicit numerical algorithm was given in [8] for each of these 
amplitudes in D = 4. 

For the practical evaluation of expectation values using numerical methods, it is often necessary to employ 
Markov chain Monte Carlo techniques [27] such as the Metropolis algorithm. For such methods it is necessary 
to find a set of moves that are ergodic — where any configuration can be transformed into any other by an 
appropriate sequence of the moves. For the present case, it can be non-trivial to find such moves due to the 
constrained nature of the dual configurations. 

To describe the ergodic moves constructed for the present work, we must first identify the variables 
present at each vertex and edge, and the nature of the constraints they obey. Each plaquette carries a 
variable whose values are quantum numbers (we shall refer to them as "spins") that label the irreducible 
unitary representations of the group G. Each edge carries an intertwiner: a G-invariant map between the 
irreducible representations assigned to the incident plaquettes. 

In the present case of a four-dimensional hypercubic lattice, we have six incident plaquettes, and hence 
the intertwiner must be an invariant map from a six-fold tensor product to the trivial representation. Due 
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to the self-duality of SU (2) irreps, this is equivalent to an invariant map between a three-fold tensor product 
of irreps to a three-fold tensor product of irreps. Intertwiners of this type form a three-dimensional vector 
space, and thus require three spin labels to resolve. The basis of intertwiners used can be represented 
diagrammatically as a splitting of the six-valent vertex into a larger network with four three- valent vertices 
and three "internal" edges. We note that there is some choice in how this splitting is carried out (equivalent 
to changing the basis of intertwiners). As we are using the vertex amplitude described in [8], we shall adopt 
the corresponding choice of splitting, which is shown below in Figure [3l The three spin labels appear as i\, 
i2, and 13 in Figure [3J The reader can readily check that these form an orthogonal basis by taking the inner 
product of two different labellings of the diagram by 11,12,13. 




Figure 3. Chosen splitting of the six-valent lattice edge and spin network vertex 

On a four-dimensional hypercubic lattice, there are 24 plaquettes incident to each vertex and 8 edges. 
This leads to a vertex amplitude that depends on 24 spin labels coming from the plaquettes and 8 x 3 = 24 
spin labels coming from the edges, in total a function of 48 spins. The approach of fixing a basis and 
explicitly carrying out a tensor contraction would be extraordinarily difficult to compute, with the cost 
scaling prohibitively with the spin [TD] . However, as mentioned above, an algorithm that scales as 0(j 4 ) has 
recently been given in [S]. 

Having described the nature of the variables, it remains to define a set of crgodic moves. To provide 
ergodicity, several distinct types of moves are applied. As with the D = 3 algorithm for SU(2) constructed 
in [JJ, these fall into three distinct classes — homology moves, edge moves, and cube moves. 

Definition A.l (Homology move). Select a plane of plaquettes that spans the entire lattice (and thus forms 
a closed surface due to the periodic boundary conditions). With equal probability, choose to decrease or 
increase a fundamental unit of spin and apply this change to every plaquette in the plane. 

Due to the periodic boundary conditions, these planes form closed surfaces that wrap the 4-torus. All 
admissible deformations of these surfaces with the same globally wrapped topology are generated by the 
plaquette and edge moves. The present algorithm uses a new approach to edge and plaquette moves that 
generally improves the acceptance ratio compared to that of [7] . Recall that given any plaquette labelling, 
the allowed intertwiner labels at each edge vary within a fixed range (non-empty if the plaquette labelling 
is admissible). Rather than storing and applying moves to the intertwiner label itself, we store the offset of 
the intertwiner relative to the minimum value in the admissible range. 

Specifically, given an edge splitting and labelling of incident plaquette spins jx,---,je, the admissible 
intertwiner labellings 23 are given by the intersection of the ranges defined by the triangle inequalities. 
For example, i\ G [|j5 — Jeljjs + je]- For each edge, a relative intertwiner label i r for each of the three 
intertwiners is stored. In computing the amplitude of the configuration, the intertwiner used is i mm + i r is 
used, where i mm is the smallest admissible intertwiner label. 

We mention in passing that for certain special configurations, the amplitude may be zero despite the fact 
that all admissibility conditions are satisfied. As was the case for the D = 3 form of this algorithm [7J, we 
assume such configurations are sufficiently isolated from each other that they do not divide the configuration 
space into disconnected regions that can't be traversed by the moves. 

Definition A. 2 (Cube move). A fundamental 3-cube of the lattice is selected at random. For each plaquette 
in the boundary of this cube, a random choice is made to either increase or decrease the spin by a single 
fundamental unit of charge. 
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Unlike the cube moves in [7J, no direct changes to the intertwiner labels is needed for compatibility. The 
intertwiners labels are implicitly changed "automatically" by virtue of the admissible ranges being modified 
by the change in plaquette label. Except when the spin is at the top of the admissible range and the plaquette 
moves reduce this range, a change in plaquette labels will implicitly give a new intertwiner labelling that 
is admissible. Another benefit of this approach is that it is actually easier to implement, as the algorithm 
doesn't have keep track of the different ways in which intertwiners sit inside cubes (as was done in [7J), but 
simply which internal edges are connected to which plaquettes. 

Definition A. 3 (Edge move). An edge e and intertwiner label i associated with that edge are randomly 
selected. The relative intertwiner labelling i r is increased or decreased by two units of fundamental spin. 

We recall the change is by two units of spin due to the parity constraint; a single unit change in either 
direction will always be inadmissible. 

In summary, the operation of the algorithm can be described as follows. An edge, cube, or homology move 
type is selected at random and the proposed move is applied. The contribution to the amplitude from those 
vertices, edges, and plaquettes effected by the move is computed before and after the move and the ratio of 
these local amplitudes is used to determine acceptance or rejection of the move by the standard Metropolis 
decision procedure. 
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